##########################################
## "Voters get what they want"          ##
##########################################
## Heinrich, Kobayashi, Long            ##
##########################################

##          
## Tables for Top40/Lowest40/Top20/Lowest20 subset models
## Script 05
## June 29, 2017


## Load estimated models
########################

for(m in 1:2)
{
  which_news_cut <- list(c(1, 3), c(2, 4))
  for(j in 1:2)
  {
    if(j == 1) numbers_to_get <- paste0("t-", as.numeric(rownames(subset(ABCs$FM, B %in% which_news_cut[[m]] & C == 1))), ".Rdata")
    if(j == 2) numbers_to_get <- paste0("t-", as.numeric(rownames(subset(ABCs$FM, B %in% which_news_cut[[m]] & C == 2))), ".Rdata")
    
    all_estimates <- list.files(path="output/analysis1_subset/", pattern="Output", full.names=TRUE)
    to_get <- c()
    for(i in 1:length(numbers_to_get)) to_get <- c(to_get, all_estimates[grepl(x=all_estimates, pattern=numbers_to_get[i], fixed=TRUE)])
    
    mods_ci <- mods_coefs <- mods_tmp <- vector("list", length(to_get))
    mods_vcv_1 <- mods_vcv_2 <- matrix("", length(to_get), 2)
    for(i in 1:length(to_get))
    {
      load(to_get[i])
      
      tmp_coef <- as.matrix(ldply(.data=output$Model, .fun=function(x) x$Sol))
      tmp_vcv <- ldply(.data=output$Model, .fun=function(x) x$VCV)
      
      mods_tmp[[i]] <- output$lmer_model
      mods_coefs[[i]] <- colMeans(tmp_coef)
      mods_ci[[i]] <- colQuantiles(tmp_coef, probs=c(.025, .975))
      mods_vcv_1[i,1] <- round(mean(tmp_vcv[,1]), di=2)
      mods_vcv_2[i,1] <- round(mean(tmp_vcv[,2]), di=2)
      mods_vcv_1[i,2] <- paste0("(", round(quantile(tmp_vcv[,1], probs=.025,), di=2),
                                ", ", round(quantile(tmp_vcv[,1], probs=.975,), di=2), ")")
      mods_vcv_2[i,2] <- paste0("(", round(quantile(tmp_vcv[,2], probs=.025,), di=2),
                                ", ", round(quantile(tmp_vcv[,2], probs=.975,), di=2), ")")
    }
    
    
    ## Long version of table
    ########################
    table_new <- stargazer(type="latex", style="ajps",
                           mods_tmp, coef=mods_coefs, ci.custom=mods_ci, 
                           ci=TRUE,
                           dep.var.labels="Total aid per capita",
                           column.labels=rep(c("High news", "Low news"), each=2),
                           label="", 
                           p.auto=FALSE, report="vcs", table.placement="t",
                           font.size="tiny", model.numbers=FALSE, model.names=FALSE,
                           covariate.labels=c("Human rights violations", "Economic hierarchy",  "Violations x economic hierarchy", 
                                              "Security hierarchy",
                                              "Violations x security hierarchy",
                                              "Lagged total aid", "Democracy", "Global aid flows",
                                              "GDP per capita", "Population", "Trade", "Alliance", "Socialist",
                                              "Cold War", "Cold War x socialist", "War", 
                                              "Ally neighbor", "UN voting similarity", ifelse(j == 1, "News (Leader)", "News (human rights)"),
                                              "Refugees", "Constant"),
                           digits=2, notes="", df=FALSE, omit.stat=c("ll", "adj.rsq", "f", "rsq"))
    
    which_row <- which(grepl(x=table_new, pattern="AIC") == TRUE)
    table_new <- c(table_new[1:(which_row-1)],  
                   paste("Residual SE & ", paste(mods_vcv_1[,1], collapse="  & "), "\\\\", sep=""),
                   paste("& ", paste(mods_vcv_1[,2], collapse="  & "), "\\\\", sep=""),
                   "  & &  \\\\ ",
                   paste("Random effect SE & ", paste(mods_vcv_2[,1], collapse="  & "), "\\\\", sep=""),
                   paste(" & ", paste(mods_vcv_2[,2], collapse="  & "), "\\\\", sep=""),
                   "\\hline \\\\[-1.8ex] ")
    table_new <- table_new[-5]
    table_new <- table_new[-5]
    
    table_new <- str_replace(string=table_new, pattern="\\[t\\] ", replacement="\\[!ht\\] ")
    if(j == 1 & m == 1)
    {    
      table_new <- c(table_new, "\\end{tabular} ",
                     "\\caption{\\textbf{Posterior summaries for main models (Leader news variable, $\\leq 40^{th}$/ $\\geq 60^{th}$ percentile cutoffs for news subsets).} Each column depicts the results for 
                        a separate model based on the respective hierarchy measure. The standalone number is the median posterior
                        estimate, the range below it the 95\\% central credible interval.}",
                     "\\end{table} ")
      
      table_new[length(table_new)-1] <- paste0(table_new[length(table_new)-1], "\\label{TabA1:SubsetNYTleaderLong}")
      write(table_new, file="output/tables/table_a1_subset_NYTleader_long.tex")
    }
    if(j == 2 & m == 1)
    {
      table_new <- c(table_new, "\\end{tabular} ",
                     "\\caption{\\textbf{Posterior summaries for main models (Human rights news variable, $\\leq 40^{th}$/ $\\geq 60^{th}$ percentile cutoffs for news subsets).} Each column depicts the results for 
                     a separate model based on the respective hierarchy measure. The standalone number is the median posterior
                     estimate, the range below it the 95\\% central credible interval.}",
                     "\\end{table} ")
      
      table_new[length(table_new)-1] <- paste0(table_new[length(table_new)-1], "\\label{TabA1:SubsetNYTHRLong}")
      write(table_new, file="output/tables/table_a1_subset_NYTHR_long.tex")
    }
    if(j == 1 & m == 2)
    {
      table_new <- c(table_new, "\\end{tabular} ",
                     "\\caption{\\textbf{Posterior summaries for main models (Leader news variable, $\\leq 20^{th}$/ $\\geq 20^{th}$ percentile cutoffs for news subsets).} Each column depicts the results for 
                     a separate model based on the respective hierarchy measure. The standalone number is the median posterior
                     estimate, the range below it the 95\\% central credible interval.}",
                     "\\end{table} ")
      
      table_new[length(table_new)-1] <- paste0(table_new[length(table_new)-1], "\\label{TabA1:SubsetNYTleaderLong-P20}")
      write(table_new, file="output/tables/table_a1_subset_NYTleader_long_p20.tex")
    }
    
    ## Short version
    ################
    which_rows <- c(which(grepl(x=table_new, pattern="Lagged total aid") == TRUE) - 1,
                    which(grepl(x=table_new, pattern="Constant") == TRUE) + 2)
    table_new <- table_new[c(1:which_rows[1], which_rows[2]:length(table_new))]
    table_new[5] <- "\\footnotesize "
    if(j == 1 & m == 1)
    {
      table_new <- c(table_new[1:34], "\\caption{\\textbf{Posterior summaries for main models; abbreviated table.} Each column depicts the results for 
               a separate model based on the respective hierarchy measure. The standalone number is the median posterior
               estimate, the range below it the 95\\% central credible interval. The full results are shown in Table \\ref{TabA1:SubsetNYTleaderLong}}\\label{TabA1:SubsetNYTleaderShort}",
                     table_new[36])
      write(table_new, file="output/tables/table_a1_subset_NYTleader_short.tex")
    }
    if(j == 2 & m == 1)
    {
      table_new <- c(table_new[1:34], "\\caption{\\textbf{Posterior summaries for main models; abbreviated table.} Each column depicts the results for 
                   a separate model based on the respective hierarchy measure. The standalone number is the median posterior
                   estimate, the range below it the 95\\% central credible interval. The full results are shown in Table \\ref{TabA1:SubsetNYTHRLong}}\\label{TabA1:SubsetNYTHRShort}",
                     table_new[36])
      write(table_new, file="output/tables/table_a1_subset_NYTHR_short.tex")
    }
  }
}




